Ciencia de Datos Geográficos aplicada a la Política Pública

Dr. Juan Pablo Carranza

Introducción al manejo de datos territoriales en R

Antes de comenzar

Para ayudar a la dinámica del taller se sugiere las y los participantes que descarguen R y R-Studio. Asímismo, se alienta a repasar algunos de los principios básicos de este lenguaje de programación.

• Para instalar R

• Para instalar R-Studio

• Nociones básicas de R

Objetivos

Los objetivos de este taller son los siguientes:

• Introducirse en el manejo de distintos tipos de datos geográficos para la representación territorial de la investigación social.

• Incorporar herramientas para la representación de datos geográficos en mapas estáticos y en mapas interactivos.

• Adquirir los conocimientos necesarios para continuar de manera autogestiva la exploración de las posibilidades del ecosistema R, aplicado a los propios problemas de investigación.

• Repasar de manera introductoria algunos elementos básicos de modelos de econometría espacial.

• Repasar algunas herramientas para la aplicación de técnicas de aprendizaje computacional a los fenómenos territoriales.

Durante el recorrido de este taller haremos un uso intensivo de la librería sf (pebesma2023?). Esta librería implementa en R un manejo eficiente y estandarizado de datos espaciales utilizando el modelo de geometrías de características simples (o simple features), un estándar internacional definido por ISO 19125 y desarrollado por el Open Geospatial Consortium (OGC). Este modelo es ampliamente utilizado en herramientas de software geográfico como GeoJSON, ArcGIS, QGIS, PostGIS, MySQL Spatial Extensions, y Microsoft SQL Server, facilitando la interoperabilidad y el análisis de datos geoespaciales en diversos entornos.

En R, el paquete sf proporciona una representación poderosa y flexible de datos espaciales vectoriales. Los objetos principales que maneja son extensiones de los clásicos data frames o tibbles, y están diseñados para facilitar la integración con otras herramientas de análisis de datos en R. Cada objeto del tipo sf incluye una columna especial denominada simple feature geometry list column, que almacena la geometría de cada observación en un formato tipo lista. Esta geometría se asocia con una o más variables adicionales (atributos) que describen cada entidad geográfica, lo que convierte cada fila en una simple feature o característica simple, compuesta tanto por sus atributos como por su estructura geométrica.

Además, sf introduce una serie de funciones que permiten realizar transformaciones geométricas, operaciones topológicas y análisis espaciales de manera eficiente, aprovechando bibliotecas subyacentes como GEOS, GDAL y PROJ. Este enfoque no solo facilita la manipulación de datos geoespaciales, sino que también asegura una alta compatibilidad con los formatos y herramientas geográficas más utilizados a nivel global.

Para comenzar, instalamos la librería y la cargamos a nuestro entorno de trabajo:

install.packages("sf")
library(sf)

Tipos de datos espaciales

Es importante distinguir entre los dos tipos principales de datos espaciales: datos vectoriales y datos raster. Ambos tipos de datos permiten representar información geográfica, pero lo hacen de manera muy diferente.

Datos vectoriales

Los datos vectoriales son el enfoque utilizado por la librería sf y se basan en representaciones geométricas precisas de objetos espaciales a través de puntos, líneas y polígonos. Estos son comúnmente llamados features o características simples. Los tipos principales de geometrías que manejan los datos vectoriales son:

Puntos

Representan ubicaciones específicas en un espacio de coordenadas. Cada punto está definido por un par de coordenadas (x, y) en el espacio 2D, o (x, y, z) en el espacio 3D.

Para comprender mejor la lógica, vamos a crear un punto en el espacio bidimensional de las variables x e y. En primer lugar, definimos un vector con las coordenadas que asignaremos a este punto:

coordenadas = c(5,3)
coordenadas

## [1] 5 3

Es decir, el punto estará ubicado en el plano carteriano en una localización definida por x = 5 e y = 3.

A continuación, utilizamos la función st_point para convertir a este vector en un objeto espacial, a partir de esas coordinadas definidas previamente:

punto <- st_point(coordenadas)
punto

## POINT (5 3)

Si graficamos este objeto, observamos lo siguiente:

plot(punto, axes = TRUE)

Si bien este punto está localizado en el plano cartesiano a partir de los valores 5 (para x) y 3 (para y), la lógica aplica cuando estos valores no representan ejes abstractos sino un sistema de coordenadas dado por la latitud y la longitud. En estos casos, como veremos más adelante, el punto estaría localizado sobre el espacio geográfico.

También podríamos contar con una colección de puntos dentro de un mismo objeto. Para ver esto, en primer lugar definimos un vector de valores para x junto con un vector de valores para y, haciendo:

px = c(5, 4, 2)
py = c(3, 5, 4)

Ahora, combinamos este conjunto de valores en un mismo data frame mediante la función cbind():

coordenadas <- cbind(px, py)
coordenadas

##      px py
## [1,]  5  3
## [2,]  4  5
## [3,]  2  4

Como se puede apreciar, el primer punto estará ubicado en la localización definida por los valores x = 5 e y = 3. El segundo en la localización dada por x = 4 e y = 5, y el tercer punto en la localización definida por el par ordenado x = 2 e y = 4.

Al igual que en el caso de un solo punto, convertimos a este data frame en un objeto espacial, pero ahora utilizando la función st_multipoint():

puntos <- st_multipoint(coordenadas)
puntos

## MULTIPOINT ((5 3), (4 5), (2 4))

Graficando este objeto:

plot(puntos, axes = TRUE)

Líneas

Las líneas son secuencias de puntos conectados que representan entidades lineales como caminos, ríos o vías de tren. Para comprender la lógica de este tipo de datos vectoriales, comencemos por definir un conjunto de coordenadas para la línea “l1”.

l1x <- c(0, 3, 5, 8, 10)
l1y <- c(0, 4, 4, 8, 10)
coordenadas_l1 <- cbind(l1x, l1y)
coordenadas_l1

##      l1x l1y
## [1,]   0   0
## [2,]   3   4
## [3,]   5   4
## [4,]   8   8
## [5,]  10  10

Es decir que la línea l1 partirá el origen (en donde x = 0 e y = 0), luego se dirigirá hacia el punto dado por x = 3 e y = 4, luego pasará por x = 5 e y = 4, y proseguirá su viaje hasta arribar a destino en la localización dada por el punto x = 10 e y = 10.

A partir de estas coordenadas, vamos a crear un objeto espacial utilizando la función st_linestring():

l1 <- st_linestring(coordenadas_l1)
plot(l1, axes = TRUE)

A continuación, podemos repetir el mismo procedimiento para crear una línea diferente en el mismo espacio cartesiano:

l2x <- c(2, 3, 12)
l2y <- c(4, 8, 5)
coordenadas_l2 <- cbind(l2x, l2y)
l2 <- st_linestring(coordenadas_l2)
plot(l2, axes = TRUE)

Finalmente, podemos combinar ambas líneas en un mismo objeto espacial, haciendo uso de la función st_multilinestring():

lineas <- st_multilinestring(list(l1, l2))
plot(lineas, axes = TRUE)

Polígonos

Los polígonos son representaciones formadas por líneas que encierran un área, como límites de países, lagos o áreas protegidas. Para comprender la lógica de este tipo de datos vectoriales, comencemos por definir un conjunto de coordenadas para un polígono imaginario en el plano cartesiano.

x <- c(5, 10, 10, 6, 5)
y <- c(1, 2, 5, 4, 1)
coordenadas <- cbind(px, py)
coordenadas

##      px py
## [1,]  5  3
## [2,]  4  5
## [3,]  2  4

La conformación del polígono parte de la lectura de esta lista de localizaciones en sentido antihorario. Es decir, la línea que define el contorno del polígono partirá, por ejemplo, del punto dado por x = 5 e y = 1. De allí se dirigirá hacia el punto localizado en x = 10 e y = 2. La línea continuará su viaje por cada uno de estos puntos que marcan los vértices del polígono, para arribar finalmente, de nuevo al punto de partida, dado por x = 5 e y = 1. Para convertir esta lista de coordenadas en un objeto espacial, utilizamos la función st_polygon():

poligono <- st_polygon(list(coordenadas))

Si graficamos este objeto, tenemos que:

plot(poligono, axes = T)

Colecciones geométricas

Combinaciones de puntos, líneas y/o polígonos en un solo objeto. Cada una de estas geometrías puede estar acompañada de atributos, como por ejemplo el nombre de una ciudad para un punto, la longitud de una ruta para una línea, o la población de un área para un polígono. Estos atributos permiten realizar análisis espaciales avanzados y crear mapas temáticos. Para armar un objeto espacial que consista en una colección de los puntos, líneas y polígonos construidos anteriormente, utilizamos la función st_geometrycollection().

coleccion <- st_geometrycollection(list(puntos, lineas, poligono))
plot(st_sfc(coleccion), axes = TRUE)

Datos raster

A diferencia de los datos vectoriales, los datos raster dividen el espacio geográfico en una cuadrícula regular de celdas o píxeles. Cada celda contiene un valor numérico que representa información continua o discreta, como la elevación de un terreno, la temperatura en una región o la clasificación del uso del suelo. Los datos raster son ideales para representar fenómenos que varían de manera continua en el espacio, como el clima, la geología o las imágenes satelitales.

A pesar de que la librería sf está centrada en datos vectoriales, es común trabajar en conjunto con datos raster, utilizando otras librerías en R como raster (Hijmans 2023) o terra (Hijmans 2024) para gestionar y analizar este tipo de datos. Las herramientas modernas permiten realizar análisis espaciales combinados, aprovechando tanto los datos vectoriales (por ejemplo, límites políticos o redes de transporte) como los datos raster (por ejemplo, modelos de elevación digital o imágenes satelitales).

Estructuralmente, los rásters consisten en matrices de valores. Para comprender mejor la lógica de este tipo de objetivos, en primer lugar, construimos una matriz de la siguiente manera:

matriz <- matrix(c(1, 2, 3, 4, 2, NA, 2, 2, 3, 3, 3, 1), ncol = 4, nrow = 3, byrow = TRUE)
matriz

##      [,1] [,2] [,3] [,4]
## [1,]    1    2    3    4
## [2,]    2   NA    2    2
## [3,]    3    3    3    1

A continuación, vamos a utilizar la librería raster para convertir esta matriz en un objeto espacial. Primero, instalamos la librería:

install.packages("raster")

Posteriormente la cargamos al entorno de trabajo, y utilizamos la función raster() para lograr el objetivo:

library(raster)

## Cargando paquete requerido: sp

r <- raster(matriz)
plot(r)

Como se observa en el gráfico, el raster contiene diferentes valores para cada celda. Los valores más bajos se representan en color gris, en tanto que lo valores más altos se representan en color verde. Las celdas representan la estructura de la matriz, en donde, por ejemplo, en la tercera fila (la de abajo) y la tercera columna (la de la derecha) el valor de la celda es igual a 1, mientras que en la primera fila (la de arriba) y la tercera columna el valor asociado es igual a 4. La intersección entre la segunda fila y la segunda columna de la matriz tiene asociado el valor NA, que se representa como una celda en blanco en el gráfico. En R, NA representa un valor faltante o no disponible (Not Available). Es la forma en que R indica que no tiene un valor para una observación en particular, lo que puede suceder por diversas razones, como datos incompletos, errores de entrada o valores omitidos.

Extrayendo datos geográficos abiertos

La comunidad de R ha generado diferentes herramientas para la extracción de datos libres y abiertos de naturaleza geográfica. Una excelente fuente de información de este tipo es Openstreetmap (OSM) (OpenStreetMap 2023), que consiste en un proyecto colaborativo que tiene como objetivo crear un mapa mundial libre y editable por cualquier persona. Se basa en la filosofía de crowdsourcing, lo que significa que los datos del mapa son generados y mantenidos por voluntarios de todo el mundo, quienes contribuyen con información sobre calles, edificios, parques, ríos, y muchos otros elementos geográficos. Para compensar la utilización gratuita de información de la cual gozaremos gracias al esfuerzo desinteresado de miles de personas que aportan información al mapa, recomendamos que las personas se creen un usuario en OSM y naveguen por el mapa, contribuyendo a sumar información de sus entornos y a elevar la calidad del mapa colaborativo. Esto puede realizarse en el www.openstreetmap.org/. El siguiente link contiene una guia paso a paso para comenzar a colaborar con el mapeo.

La librería osmdata (Padgham et al. 2017) permite extraer datos desde la API de OSM, por lo tanto, comenzamos por instalarla. También instalaremos la librería nominatimlite (Hernangómez 2024), que necesitaremos para definir entorno geográficos conocidos sólo a partir de su nombre:

install.packages("osmdata")
install.packages("nominatimlite")

A continuación, vamos a definir un entorno geográfico (un polígono) sobre el cuál se acotará la extracción de información de OSM. Este polígono vendrá definido por los límites administrativos de la Municipalidad de Córdoba:

library(nominatimlite)
library(osmdata)
area_de_estudio <- geo_lite_sf(address = "Cordoba, Argentina", points_only = F)

Veamos un poco más de cerca el polígono area_de_estudio:

area_de_estudio

## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -64.3101 ymin: -31.52513 xmax: -64.05725 ymax: -31.30754
## Geodetic CRS:  WGS 84
## # A tibble: 1 × 3
##   query              address                                            geometry
## * <chr>              <chr>                                         <POLYGON [°]>
## 1 Cordoba, Argentina Córdoba, Municipio de Córdoba, P… ((-64.3101 -31.36306, -6…

Ya entraremos en más detalles respecto a este objeto, pero por lo pronto, notemos que tiene la estructura de un data frame, con una sola fila en donde se observan las siguientes columnas: query, address, geometry. La última columna es la que nos interesa en este momento, dado que es una lista de coordenadas geográficas que nos permitirá visualizar este objeto en un mapa.

Para ello, vamos a hacer uso de la librería mapview (Appelhans et al. 2023), que permite la construcción de mapas interactivos de una manera muy sencilla. Instalamos esta librería:

install.packages("mapview")

Y a continuación, construimos el mapa:

library(mapview)
mapview(area_de_estudio)

Al hacer click en el polígono notamos que se abre una ventana con el resto de la información contenida en las columnas restantes. Esta es la principal ventaja de la librería sf, que permite mantener la lógica de data frames al trabajar con objetos espaciales.

Ahora que hemos acotado el área sobre la cuál nos interesa extraer información, utilizaremos la librería osmdata para hacernos de la información necesaria para ejemplificar los diferentes tipos de objetos espaciales que hemos estado estudiando.

El siguiente bloque de código crea el objeto espacial “comercios” mediante la siguiente lógica: en primer lugar se utiliza la función opq() para acceder al servicio de descarga de información de OSM (una API) definiendo el área sobre la cuál se necesitará extraer la información. Esta área la hemos definido con el polígono “area_de_estudio” que representa los límites administrativos de la ciudad de Córdoba. En segundo lugar, utilizamos la función add_osm_feature() para definir qué tipo de datos queremos extraer de OSM. Openstreetmap tiene una manera de catalogar la información que la comunidad carga, definiendo diferentes códigos para diferentes tipos de objetos. Estos códigos pueden tener un primer agrupamiento, dado por el valor “key”, y un segundo agrupamiento dado por el valor “value”. En el ejemplo, nos interesa descargar todos los comercios (shops) en el área de estudio, pero podríamos estar interesados en descarga sólo las panaderías (bakery). En este hipotético caso, deberíamos reemplazar la línea correspondiente por add_osm_feature(key = “shop”, value = “bakery”). Una exploración detallada de los diferentes criterios de catalogación de la información en OSM se puede realizar en el siguiente link.

comercios <-  opq(st_bbox(area_de_estudio)) %>%
  add_osm_feature(key = "shop") %>%
  osmdata_sf()

Veamos un poco más de cerca el objeto recientemente creado, llamado “comercios”:

comercios

## Object of class 'osmdata' with:
##                  $bbox : -31.5251273,-64.3100961,-31.307545,-64.05725
##         $overpass_call : The call submitted to the overpass API
##                  $meta : metadata including timestamp and version numbers
##            $osm_points : 'sf' Simple Features Collection with 3415 points
##             $osm_lines : NULL
##          $osm_polygons : 'sf' Simple Features Collection with 374 polygons
##        $osm_multilines : NULL
##     $osm_multipolygons : NULL

De lo anterior se desprende que este objeto contiene 3415 puntos, y 374. Esto se debe a que los usuarios pueden marcar la localización de un comercio con un punto, pero también dibujar un polígono que describa grandes superficies comerciales como los shoppings o ferias. Ya veremos cómo trabajar sobre esta dispersión de criterios, pero de momento, vamos a quedarnos sólo con los 3415 puntos.

comercios = comercios$osm_points

Veamos un poco en detalle el resultado. La función dim() nos informa sobre la dimensión del data frame resultante.

dim(comercios)

## [1] 3415  111

Como vemos, el objeto “comercios” tiene 3415 filas y 111 columnas. Es mucha información. Si utilizamos la función names() podemos ver los nombres de cada una de estas columnas:

names(comercios)

##   [1] "osm_id"                         "name"                          
##   [3] "access"                         "addr:city"                     
##   [5] "addr:country"                   "addr:floor"                    
##   [7] "addr:housenumber"               "addr:postcode"                 
##   [9] "addr:street"                    "air_conditioning"              
##  [11] "alt_name"                       "amenity"                       
##  [13] "barrier"                        "branch"                        
##  [15] "brand"                          "brand:en"                      
##  [17] "brand:wikidata"                 "brand:wikipedia"               
##  [19] "building"                       "bus"                           
##  [21] "butcher"                        "check_date"                    
##  [23] "clothes"                        "contact:email"                 
##  [25] "contact:facebook"               "contact:instagram"             
##  [27] "contact:phone"                  "contact:website"               
##  [29] "currency:ARS"                   "currency:XBT"                  
##  [31] "delivery"                       "description"                   
##  [33] "designation"                    "diet:vegetarian"               
##  [35] "drink:coffee:togo"              "drink:wine"                    
##  [37] "email"                          "entrance"                      
##  [39] "female"                         "fixme"                         
##  [41] "healthcare"                     "highway"                       
##  [43] "instagram"                      "internet_access"               
##  [45] "internet_access:fee"            "is_in"                         
##  [47] "laundry_service"                "level"                         
##  [49] "male"                           "mobile"                        
##  [51] "name:en"                        "name:es"                       
##  [53] "network"                        "not:brand:wikidata"            
##  [55] "note"                           "official_name"                 
##  [57] "opening_hours"                  "operator"                      
##  [59] "organic"                        "outdoor_seating"               
##  [61] "payment:american_express"       "payment:app"                   
##  [63] "payment:argencard"              "payment:bank_deposits"         
##  [65] "payment:bank_transfer"          "payment:cabal"                 
##  [67] "payment:cash"                   "payment:cheque"                
##  [69] "payment:contactless"            "payment:credit_cards"          
##  [71] "payment:debit_cards"            "payment:diners_club"           
##  [73] "payment:electronic_purses"      "payment:gift_card"             
##  [75] "payment:italcred"               "payment:kadicard"              
##  [77] "payment:maestro"                "payment:mastercard"            
##  [79] "payment:mastercard_contactless" "payment:mastercard_debit"      
##  [81] "payment:meal_voucher"           "payment:naranja"               
##  [83] "payment:nativa"                 "payment:provencred"            
##  [85] "payment:tarjeta_shopping"       "payment:visa"                  
##  [87] "payment:visa_debit"             "payment:visa_electron"         
##  [89] "payment:wire_transfer"          "phone"                         
##  [91] "phone:AR"                       "prepayment:cencosud"           
##  [93] "public_transport"               "ref"                           
##  [95] "ref_code"                       "repair"                        
##  [97] "second_hand"                    "self_service"                  
##  [99] "service"                        "shop"                          
## [101] "shop_1"                         "shop_2"                        
## [103] "source"                         "start_date"                    
## [105] "takeaway"                       "toilets"                       
## [107] "trade"                          "website"                       
## [109] "wheelchair"                     "wholesale"                     
## [111] "geometry"

El objeto tiene una columna de geometría (la última) que define las coordenadas geográficas de cada punto, pero también muchas columnas con información adicional, que el usuario que realiza la carga de la información puede o no completar.

Para reducir el ruido de la información resultante, vamos a quedarnos sólo con la columna “shop”, que indica el rubro de negocio al cual se dedica cada uno de los comercios:

comercios = comercios[, c("shop")]
head(comercios)

## Simple feature collection with 6 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -64.21934 ymin: -31.43437 xmax: -64.14273 ymax: -31.36709
## Geodetic CRS:  WGS 84
##                  shop                    geometry
## 643269288        mall POINT (-64.21934 -31.36709)
## 655193217      bakery POINT (-64.18184 -31.42471)
## 663243931 supermarket  POINT (-64.1477 -31.42402)
## 664087358 supermarket POINT (-64.14273 -31.43437)
## 705466795     chemist POINT (-64.18341 -31.42552)
## 716005842     alcohol POINT (-64.18664 -31.42056)

Y a continuación podemos crear un mapa interativo con los datos resultantes:

mapview(comercios, col.regions = "blue", legend = F)

Siguiendo el procedimiento anterior, podemos también extraer otros tipos de objetos espaciales. Por ejemplo, las líneas que representan calles o avenidas catalogadas como “vias secundarias”. Para ello, en la función add_osm_feature() utilizamos la key “highway” y el value “secondary”.

vias_secundarias <- opq(st_bbox(area_de_estudio)) %>%
  add_osm_feature(key = "highway", value = "secondary") %>%
  osmdata_sf()
vias_secundarias = vias_secundarias$osm_lines[, c("name")]
mapview(vias_secundarias, color = "blue", legend = F)

## Warning in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy, :
## GDAL Error 1: PROJ: proj_as_wkt: DatumEnsemble can only be exported to
## WKT2:2019

Por supuesto, también podemos extraer polígonos. En este caso, por ejemplo, extraeremos los polígonos de las escuelas localizadas dentro del área de estudio haciendo uso de la key “amenity” y el value “school”. Entonces:

escuelas <- opq(st_bbox(area_de_estudio)) %>%
  add_osm_feature(key = "amenity", value = "school") %>%
  osmdata_sf()
escuelas = escuelas$osm_polygons[, c("name")]
mapview(escuelas, color = "black", col.regions = "red", legend = F)

A continuación, exploraremos algunas estrategias para utilizar datos en formato ráster. Para ello, necesitaremos instalar previamente las librerías terra y geodata (Hijmans et al. 2024).

install.packages("terra")
install.packages("geodata")

A continuación, utilizaremos la librería geodata para descargar un modelo digital de elevación en formato ráster. Esta librería permite acceder a informacion libre generada en el marco de la mision Shuttle Radar Topography Mission (SRTM), que determina la elevacion del terreno (en metros sobre el nivel de mar) a traves de un radar incorporado a un satelite. Dado que la informacion se encuentra “precargada” en mozaicos que cubren diferentes extensiones geográficas, se le debe indicar a la función elevation_3s() las coordenadas geográficas de un punto, para que a partir de esta información se puede determinar qué mosaico se va a descargar.

Ahora bien, hasta aquí hemos definido el área de estudio como un polígono definido por el objeto area_de_estudio. Como ya sabemos, este polígono no tiene sólo una, sino varias coordenadas que, leidas en sentido antihorario, definen los límites del área. Y la función requiere de sólo un par de coordenadas, es decir, de un punto. Para sortear este inconvneiente, utilizaremos la función st_centroid() de la librería sf, que permite calcular el centroide de un polígono, entendido como el punto central de un polígono que se calcula a partir del promedio de las coordinadas geográficas que lo definen. Por lo tanto:

coordenadas = st_centroid(area_de_estudio)

A continuación, a este nuevo objeto espacial que representa el centroide del polígono del área de estudio, le extraeremos las coordenadas mediante la función de sf st_coordinates:

coordenadas = st_coordinates(coordenadas)
coordenadas

##              X        Y
## [1,] -64.18378 -31.4166

En donde X representa la longitud (este-oeste) e Y representa la latitud (norte-sur). Ahora sí, estamos en condiciones de utilizar estas dos coordenadas geográficas para indicarle a la función elevation_3s() de la librería geodata el lugar para el cual queremos descargar la información sobre la elevación del terreno. La función tiene varios arguementos: el primero (lon) requiere que se informe la longitud calculada recién (el valor de X), y para asgnar ese valor tomamos el primer elemento del objeto coordenadas; para informar la latitud (lat) hacemos lo mismo pero con el segundo elemento del objeto coordenadas; finalmente, el argumento “path” requiere informar la ruta en la cual se guardará el archivo a descargar, en donde indicamos que se almacene en un directorio temporal (esto lo pueden modificar ustedes, detallando un directorio específico en sus computadoras):

library(geodata)

## Cargando paquete requerido: terra

## terra 1.8.9

library(terra)
elevacion <- elevation_3s(lon = coordenadas[[1]], lat = coordenadas[[2]], path=tempdir() )

Podemos graficar el ráster descargado:

plot(elevacion)

También podemos hacer un mapa interactivo utilizando la librería mapview, pero aparece un pequeño inconveniente. El objeto “elevación” es de tipo SpatRaster, que es el formato con el que la librería terra trabaja los rásters. Dado que mapview es una librería diseñada con anterioridad a la aparición de la librería terra, la interacción de ambas no es óptima. Para sortear este problema, convertimos el objeto “elevación” de SpatRaster a un ráster clásico de la librería raster.

mapview(raster(elevacion))

Sin embargo, el ráster descargado es mucho más amplio que el área de estudio definida. Para adecuar esta información hacemos uso de la función crop() de la librería terra, de la siguiente manera:

cordoba = crop(elevacion, area_de_estudio)
plot(cordoba)

Y al igual que recién, podemos hacer un mapa interactivo con la salvedad de convertir el objeto “cordoba” al formato correpondiente:

mapview(raster(cordoba))

Cargar datos de manera local (desde la PC)

Si bien hasta ahora hemos descargado una amplia gama de archivos con datos espaciales directamente desde internet, hay muchas fuentes de información que no tienen sevicios para conectarse directamente y requieren de la descarga previa de algún archivo. Los archivos en formato geográfico más conmunes son aquellos con las extensiones gkpg, kmz, geojson o shp. Sin embargo, cualquier archivo que tenga una columna con la lista de coordenadas geográficas de cada fila puede convertirse en un objeto espacial, incluso una simple hoja de cálculo.

Vamos a utilizar la función st_read() de la librería sf para cargar al entorno de trabajo el archivo barrios_cordoba.gpkg, creando un nuevo objeto llamado “barrios”, y generamos el mapa interactivo correspondiente:

barrios <- st_read("barrios_cordoba.gpkg")

## Reading layer `barrios_cordoba' from data source 
##   `/home/juan/Documentos/GitHub/spatial_data_science/barrios_cordoba.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 485 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XYZ
## Bounding box:  xmin: -64.30993 ymin: -31.52498 xmax: -64.05738 ymax: -31.30852
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

mapview(barrios)

Como se mencionó, existen muchos portales que ofrecen descargas de información geográfica que puede ser luego cargada de esta manera al entorno de trabajo en R. Recomendamos revisar los portales poblaciones.org, datos abiertos de la Municipalidad de Córdoba, el portal de la Infraestructura de Datos Espaciales de la Provincia de Córdoba, el portal del INDEC. Y hay muchos otros.

Operaciones con datos espaciales

A continuación, analizaremos algunas funciones de la librería sf que permiten realizar diferentes operaciones con datos espaciales. Comenzaremos por la realización de uniones entre diferentes objetos espaciales. La lógica de esta operación implica, por ejemplo, identificar los atributos de objetos que comparten una localización o localizaciones próximas.Supongamos que nos interesa imputar a cada escuela el barrio al que pertenece. Para ello, utilizaremos la función st_join(). Esta función tiene 3 argumentos. En primer lugar se indica el objeto espacial sobre el que se desea agregar un atributo correspondiente a otro objeto espacial, en este caso, las escuelas. En segundo lugar se indica el objeto espacial desde dónde se desea transferir el atributo correspondiente, en este caso, el nombre del barrio. En tercer lugar, se debe indicar el criterio de union. En este caso, se indica la opción “st_intersects” que establece que ambos objetos deben unirse siempre que se intersecten en algún punto. Hay un cuarto argumento, que no es obligatorio utilizar, pero suele ser útil: se trata de a opción “largest”. En el caso de asignarle el valor “TRUE”, y en el frecuente caso de que una feature de un objeto se intersecte con dos o más features del otro objeto (por ejemplo, que una escuela esté ubicada en el límite de dos barrios, con una parte en uno y otra parte en otro) esta opción vinculará los datos de aquellas que compartan la mayor proporción del espacio.

escuelas = st_join(escuelas, barrios[,c("description")], join = st_intersects, largest = TRUE)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

head(escuelas)

## Simple feature collection with 6 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -64.27902 ymin: -31.45692 xmax: -64.1732 ymax: -31.36073
## Geodetic CRS:  WGS 84
##                                                   name
## 67868340                      Instituto Nuestra Señora
## 83552775           Instituto Técnico Salesiano Villada
## 143976290 Escuela Superior de Comercio Manuel Belgrano
## 154249743          Instituto Técnico Salesiano Villada
## 154644719                     Colegio San Buenaventura
## 184007376                    IPEM 8 Manuel Reyes Reina
##                          description                       geometry
## 67868340  AMPLIACION JARDIN ESPINOSA POLYGON ((-64.17479 -31.454...
## 83552775             VALLE ESCONDIDO POLYGON ((-64.27749 -31.363...
## 143976290                    ALBERDI POLYGON ((-64.20047 -31.405...
## 154249743            VALLE ESCONDIDO POLYGON ((-64.27666 -31.362...
## 154644719         ALTOS DE SANTA ANA POLYGON ((-64.25112 -31.402...
## 184007376                 ALTO VERDE POLYGON ((-64.21076 -31.369...

Como vemos, se ha agregado una columna que informa el barrio al que pertenece cada escuela.

Ahora, con esta nueva información, podemos realizar algunas agregaciones no espaciales, pero que tendrán correlato a nivel geográfico. Por ejemplo, supongamos que nos interesa saber la cantidad de escuelas que hay en cada barrio. Para ello, recuperando lo estudiado en módulos anteriores, hacemos uso de las funciones group_by() y summarise() de la librería dplyr (Wickham et al. 2023). Antes de proceder a calcular ese resumen de cantidad de escuelas por barrio, vamos a remover los atributos espaciales del objeto “escuelas”. Esto lo haremos utilizando la función st_drop_geometry(), que básicamente convierte a data frame todo objeto espacial. Entonces

library(dplyr)
tabla <- st_drop_geometry(escuelas) %>% 
  group_by(description) %>% 
  summarise(cantidad = n()) %>% 
  arrange(desc(cantidad))
head(tabla)

## # A tibble: 6 × 2
##   description         cantidad
##   <chr>                  <int>
## 1 ALBERDI                    9
## 2 CERRO DE LAS ROSAS         7
## 3 GENERAL PAZ                7
## 4 ARGUELLO                   6
## 5 ARGUELLO NORTE             6
## 6 VILLA AZALAIZ OESTE        6

Ahora supongamos que nos interesa confeccionar un mapa que informe sobre la cantidad de escuelas en cada barrio. Vamos a hacer una unión no espacial entre este objeto no espacial llamado “tabla” con el objeto espacial llamado “barrios”. La llave para unir ambas tablas será el campo description que se encuentra presente en ambas. Por lo tanto:

barrios = left_join(barrios, tabla)

## Joining with `by = join_by(description)`

Ahora sí, confeccionamos un mapa interactivo que informe sobre la cantidad de escuelas en cada barrio de la ciudad. Como es usual, hacemos uso de la función mapview(), agregando ahora el argumento “zcol” que indica sobre qué columna se pretende agregar una paleta de colores a la representación gráfica.

mapview(barrios, zcol = "cantidad")

Es importante notar que muchos barrios (330 de un total de 485) tienen NA en el campo “cantidad”. Esto se debe a que, del total de escuelas descargadas desde OSM, ninguna se encuentra localizada dentro del polígono de cada uno de estos barrios. Conviene recordar en este punto que la fuente de información desde la cual se extrajo el dato de escuelas no es oficial, sino que se trata de un mapa colaborativo en el que, seguramente, faltan varias por agregar (aprovachamos para reforzar la necesidad de colaborar con OSM para mejorar la calidad de esta información).

Sobre el mapa anterior podemos agregar las escuelas, para asegurarnos de que, efectivamente, no hay ninguna en los barrios que aparecen con valores NA (en gris en el mapa anterior):

mapview(barrios, zcol = "cantidad") +
  mapview(escuelas, col.regions = "red")

Otra operación espacial interesante consiste en la unión de diferentes features. Por ejemplo, supongamos que nos interesa unir los barrios Centro, Nueva Córdoba y Alberdi. Lo primero que vamos a hacer es filtrar esos barrios utilizando la función filter().

barrios <- barrios %>% filter(description %in% c("CENTRO","NUEVA CORDOBA","ALBERDI"))
mapview(barrios)

Como vemos, el objeto espacial “barrios” tiene ahora sólo 3 filas o features, una para cada barrio. Vamos a unir estas 3 features en un solo polígono, utilizando la función st_union() de la librería sf.

barrios = st_union(barrios)
mapview(barrios)

Como vemos, ahora el objeto tiene sólo un polígono, resultante de la union de los 3 barrios anteriores.

Supongamos que nos interesa ver sólo las escuelas que se encuentran dentro de este nuevo polígono. Para ello, utilizaremos la función st_intersection() de la librería sf. Esta función tiene sólo dos argumentos: el primero consiste en la capa sobre la que se quiere seleccionar las features (las escuelas en nuestro ejemplo), y el segundo representa el área sobre la cuál se desea realizar la selección (el polígono de los 3 barrios unidos). Entonces:

escuelas = st_intersection(escuelas, barrios)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

mapview(escuelas)

El objeto espacial resultante contiene ahora exclusivamente las 15 que se encuentran localizadas dentro de los 3 barrios seleccionados.

Introducción a la econometría espacial

Los modelos de regresión lineal tradicionales son inadecuados cuando los supuestos para su implementación dejan de ser validos. Si las unidades espaciales están caracterizadas por relaciones de dependencia o autocorrelación espacial, estas deben ser tenidas en cuenta para evitar problemas en la especificación del modelo.

(Anselin 1988) plantea la distinción entre autocorrelación espacial o dependencia espacial y heterogeneidad espacial o no estacionariedad espacial. Esta última plantea la posibilidad de que los parámetros sean variables en el espacio. Se focalizará en la dependencia o autocorrelación espacial, la cual existe cuando se pueden apreciar similitudes entre observaciones en una misma localización, es decir cuando el valor de la observación tiene una estrecha relación con la de sus vecinos (ley de Tobler). Considerando el foco en la dependencia espacial, se verá a continuación cómo se puede “atrapar” esta dependencia con modelos de regresión espacial.

Los modelos espaciales pueden ser clasificados en función de los tipos de interacciones espaciales que poseen. Estas pueden ser:

• Interacciones endógenas, donde la decisión económica de un agente o zona geográfica va a depender de la decisión de sus vecinos.

• Interacciones exógenas, donde la decisión económica de un agente o zona geográfica va a depender de las características observables de sus vecinos.

• Correlación espacial como concecuencia de caracteristicas no observadas.

En la siguiente ecuación se escribe el modelo de (Manski 1993) escrito en forma matricial:

\[Y = \rho W Y + \beta X + \theta W X + \mu\] \[ \mu = \lambda W \mu + \epsilon\] Donde \(\beta\) son los parámetros (a estimar) de las variables exógenas explicativas, \(\rho\) es el el efecto de interacción endógena (efecto autorregresivo espacial), \(\rho\) son los efectos de interacción exógena y \(\lambda\) son los efectos de correlación espacial de errores (autocorrelación espacial).

La matriz \(W_{ij}\) juega un rol central en el análisis. Se conoce como matriz de vecindario y tiene una dimensión igual a \((n\) x \(n)\) siendo \(n\) la cantidad de observaciones. En su versión más simple cada escalar asume el valor 1 si la observación \(i\) es vecina de la observación \(j\), y cero en otro otro caso. La diagonal princial de la matriz es, por definición, igual a cero (dado que una observación no puede ser vecina de sí misma). La especificación de la matriz se puede hacer más compleja, por ejemplo, dividiendo cada elemento por la distancia entre pares ordenados de observaciones \((i, j)\), de manera tal de otorgar un mayor peso a las observaciones más próximas. Existen múltiples criterios para definir la “vecindad” entre observaciones. Uno de los criterios más utilizados, y el más simple, consiste en definir un vecindario a partir de un radio de \(m\) metros de cada observación.

Resulta conveniente resaltar que resolviendo algebraicamente, por ejemplo, \(\rho W Y\), se obtiene un vector que indicar el promedio (ponderado si se incorporó la distancia a la especificación de la matriz \(W\)) del valor de \(Y\) en el vecindario.

El modelo no puede ser identificable de esta forma (\(\rho\),\(\lambda\),\(\beta\),\(\theta\)) no pueden ser estimados al mismo de manera simultánea dados los grados de libertad). Para que este modelo sea identificable se puede elegir entre dos opciones: utilizar matrices de vecindario distintas o eliminar alguna de las tres formas de correlación espacial enunciadas. La última opción es la preferida empíricamnte.

Por lo tanto, asumiendo diferentes opciones para incorporar la dependencia espacial al modelo, la literatura provee de las siguientes especificaciones:

• SLM (Standard Linear Regression Model): Es el modelo lineal sin dependencia espacial, es decir que: \(\rho=0\), \(\theta=0\) y \(\lambda=0\):

\[Y = \beta X + \epsilon\]

\[\epsilon \sim N(0, \sigma^2 I_n)\]

• SEM (Spatial Error Model): Si tanto: \(\rho=0\) como \(\theta=0\) se considera la dependencia espacial sólo en el término de error:

\[Y = \beta X + \mu\]

\[\mu = \lambda W \mu + \epsilon\]

• SAR (Spatial Autoregressive Model): Si tanto: \(\lambda=0\) como \(\theta=0\) se considera la dependencia espacial sólo en la variable dependiente:

\[Y = \rho W Y + \beta X + \epsilon\]

\[\epsilon \sim N(0, \sigma^2 I_n)\]

• SLX (Spatially Lagged X): Si tanto: \(\lambda=0\) como \(\rho=0\) se considera la dependencia espacial sólo en las variables independientes:

\[Y = \beta X + \theta W X + \epsilon\]

\[\epsilon \sim N(0, \sigma^2 I_n)\]

• SAC (Spatially Autoregressive Confused): Si sólo \(\theta=0\) se considera la dependencia espacial tanto en la variable dependiente como en el término de error:

\[Y = \rho W Y + \beta X + \mu\]

\[\mu = \lambda W \mu + \epsilon\]

• SDM (Spatial Durbin Model): Si sólo \(\lambda=0\) se considera autocorrelación espacial en las variables independientes y en la variable dependiente.

\[Y = \rho W Y + \beta X + \theta W X + \epsilon\]

\[\epsilon \sim N(0, \sigma^2 I_n)\]

• SDEM (Spatial Durbin Error Model): Si sólo \(\rho=0\) e asume que no hay interacción endógena y que el énfasis se ubica en las externalidades de los vecinos y el error:

\[Y = \beta X + \theta W X + \mu\]

\[\mu = \lambda W \mu + \epsilon\]

Selección de Modelos

Para elegir qué modelo aplicar existen aproximaciones prácticas basadas en los supuestos de que existe una matriz de vecindario conocida, que las variables explicativas son exógenas y que \(\epsilon \sim N(0, \sigma^2 I_n)\).

Sin embargo, es muy importante el marco teórico del problema para realizar la elección de modelos. Las metodologías que se presentan a continuación sirven como guías para la elección del modelo, pero siempre debe haber una revisión bibliográfica que sustente estas relaciones.

Para identificar los modelos se presentan tres posibles procedimientos a seguis: un enfoque bottom-up (de abajo hacia arriba), un enfoque top-down (de arriba hacia abajo), enfoque combiando y por último el enfoque de reducción por modelos anidados.

Multiplicadores de Lagrange

Los multiplicadores de lagrange (Anselin 1995) en este contexto se utilizan para probar la hipótesis nula de que la dependencia espacial no existe en los residuos o en la variable dependiente. Es decir, se pueden plantear dos test:

• SEM vs SLM: donde H0: \(\lambda=0\)
• SAR vs SLM: donde H0: \(\rho=0\)

Es decir que estos multiplicadores prueban si el modelo mejora su ajuste si se utiliza un ajuste por dependencia espacial en lugar de un modelo linal. Para concluir que es mejor utilizar un modelo espacial de los tipos mencionados se debe rechazar la hipótesis nula planteada. Es decir, el test debe mostrar un p-valor < \(\alpha\) .

Se pueden utilizar (de hecho se recomienda hacerlo) los multiplicadores de lagrange robustos, ya que permiten filtrar falsos positivos, por lo que las conclusiones serán más robustas. La iterpretación es análoga a la realizada anteriormente.

Bottom-up Approach (abajo-arriba)

Este enfoque consiste en comenzar con el modelo lineal no espacial y utilizar el test de multiplicadores de Lagrange para elegir entre el modelo SAR, SEM o el no espacial. Adicionalmente, si se encuentra que los parámetros espaciales (\(\rho\) o \(\lambda\)) son ambos significativamente distintos de cero, se puede optar un un modelo SAC. Como se puede apreciar, en este enfoque no se considera la posibilidad de elegir entre modelos que incluyan la dependencia espacial en las variables independientes.

Top-Down Approach (arriba-abajo)

Para este enfoque se comienza desde un modelo SDM (dependencia espacial exógena y endógena) y se va deduciendo el modelo más adecuado.

Abordaje combinado

Como el título lo indica, es una combinación de los enfoques ateriores. Se comienza con la estimación de un modelo lineal (bottom-up) y ante la presencia de interacciones espaciales (\(\lambda \neq 0\) o $) se tiene en cuenta la posibilidad de la dependencia exógena. Hay que considerar que si se se decide estimar un modelo SAC, ya no queda posibilidad de seguir estudiando una dependencia espacial adicional en las variables independientes porque se pierde la identificabilidad del modelo.

¿Cómo hacerlo en R?

Suponga que cuenta con una base de datos con el precio por metro cuadrado de la tierra urbana, junto con algunas características que informan sobre la localización de cada lote respecto de diferentes hitos urbanos y algunas características del entorno. Las variables son las siguientes:

• valor = precio por metro cuadrado de la tierra urbana, en dólares estadounidenses.
• dist_bancos: distancia al banco más cercano.
• dist_barrio_cerrado: distancia al barrio cerrado más cercano.
• dist_barrios_populares: distancia al barrio popular o asentamiento informal más cercano.
• dist_industrias: distancia a la industria más cercana.
• dist_vias_primarias: distancia a la vía primaria más cercana.
• dist_vias_secundarias: distancia a la vía secundaria más cercana.
• entorno_cuentas: cantidad de unidades habitacionales en un entorno de 500 metros.
• entorno_edificaciones: porcentaje del territorio que se encuentra edificado en un entorno de 500 metros.
• entorno_fot: promedio del cociente entre la cantidad de metros cuadrados edificados en cada parcela y su superficie, en un entorno de 500 metros.
• entorno_parcelas: cantidad de parcelas en un entorno de 500 metros.
• entorno_plazas: porcentaje del territorio que se encuentra ocupado por plazas o parques en un entorno de 500 metros.

load("datos.rda")
summary(datos)

##      valor         dist_bancos       dist_barrio_cerrado dist_barrios_populares
##  Min.   :   6.0   Min.   :   11.22   Min.   :    0.0     Min.   :    0.0       
##  1st Qu.:  33.0   1st Qu.: 1126.47   1st Qu.:  633.3     1st Qu.:  913.9       
##  Median :  66.0   Median : 2194.46   Median : 2076.8     Median : 1696.0       
##  Mean   : 150.9   Mean   : 3133.55   Mean   : 3877.8     Mean   : 3276.5       
##  3rd Qu.: 133.0   3rd Qu.: 3639.13   3rd Qu.: 4999.1     3rd Qu.: 3412.8       
##  Max.   :5500.0   Max.   :19385.13   Max.   :19631.1     Max.   :19453.5       
##  dist_industrias   dist_vias_primarias dist_vias_secundarias entorno_cuentas  
##  Min.   :    0.0   Min.   :   13.24    Min.   :   7.264      Min.   :    1.0  
##  1st Qu.:  912.3   1st Qu.:  422.85    1st Qu.: 183.888      1st Qu.:  331.0  
##  Median : 1637.2   Median : 1103.89    Median : 470.689      Median :  605.0  
##  Mean   : 2321.5   Mean   : 1889.95    Mean   : 811.212      Mean   :  940.3  
##  3rd Qu.: 2934.7   3rd Qu.: 2583.23    3rd Qu.: 994.556      3rd Qu.:  974.5  
##  Max.   :12429.0   Max.   :11969.14    Max.   :6585.321      Max.   :25206.0  
##  entorno_edificaciones  entorno_fot       entorno_parcelas entorno_plazas    
##  Min.   : 0.00         Min.   :0.000061   Min.   :   1.0   Min.   :0.000000  
##  1st Qu.:10.59         1st Qu.:0.061616   1st Qu.: 311.5   1st Qu.:0.000000  
##  Median :21.06         Median :0.123495   Median : 565.0   Median :0.006277  
##  Mean   :22.56         Mean   :0.212486   Mean   : 620.3   Mean   :0.014353  
##  3rd Qu.:35.27         3rd Qu.:0.266695   3rd Qu.: 846.0   3rd Qu.:0.021351  
##  Max.   :48.78         Max.   :3.652536   Max.   :1667.0   Max.   :0.229007  
##             geom     
##  POINT        :1139  
##  epsg:22174   :   0  
##  +proj=tmer...:   0  
##                      
##                      
## 

Si interesa realizar un análisis de regresión para estudiar el efecto de cada una de estas variables independientes sobre el precio por metro cuadrado de la tierra urbana, resulta conveniente tomar logaritmos naturales tanto en la variable dependiente como en las independientes, para estimar elasticidades y semi-elasticidades, de manera tal de facilitar la comprensión de los impactos.

Por ejemplo, suponga que la variable dependiente \(Y\) y la variable independiente \(X\) son continuas. En tal caso, al hacer:

\[ln(Y) = \beta \ ln(X)\] Y tomando el diferencial total de la expresión anterior:

\[\frac{1}{Y} \delta Y = \beta \frac{1}{X} \delta X\] Reordenando, el parámetro \(\beta\) es la elasticidad de un cambio en \(X\) sobre \(Y\).

De igual manera, si resutla que \(X\) es una variable dicotómica que asume los valores 0 o 1 (variable dummy), la diferencial calculada anteriormente resulta:

\[\frac{1}{Y} \delta Y = \beta \delta X\] En donde \(\beta\) es la semi-elasticidad que informa el efecto en \(Y\) de la ocurrencia observada en \(X\).

Por lo tanto, teniendo en cuenta lo anterior, se estima un modelo conforme a la fórmula siguiente:

form = log(valor) ~ log(dist_bancos) + log(dist_barrio_cerrado) + log(dist_barrios_populares) + log(dist_industrias) + log(dist_vias_primarias) + log(dist_vias_secundarias) + entorno_cuentas + log(entorno_edificaciones) + log(entorno_fot) + entorno_parcelas + log(entorno_plazas)

Se estima el modelo por mínimos cuadrados ordinarios:

ols = lm(form, datos)
summary(ols)

## 
## Call:
## lm(formula = form, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.2499 -0.3347  0.0057  0.2917  3.9967 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  8.863e+00  1.973e-01  44.931  < 2e-16 ***
## log(dist_bancos)            -1.786e-01  2.011e-02  -8.881  < 2e-16 ***
## log(dist_barrio_cerrado)    -2.933e-02  1.854e-03 -15.824  < 2e-16 ***
## log(dist_barrios_populares)  9.267e-03  9.755e-03   0.950   0.3423    
## log(dist_industrias)        -1.286e-02  6.245e-03  -2.059   0.0397 *  
## log(dist_vias_primarias)    -1.040e-01  1.279e-02  -8.133 1.09e-15 ***
## log(dist_vias_secundarias)  -1.015e-01  1.392e-02  -7.292 5.76e-13 ***
## entorno_cuentas              1.375e-04  9.427e-06  14.583  < 2e-16 ***
## log(entorno_edificaciones)  -1.182e-01  7.541e-03 -15.670  < 2e-16 ***
## log(entorno_fot)             5.909e-01  2.568e-02  23.015  < 2e-16 ***
## entorno_parcelas            -5.791e-04  5.701e-05 -10.158  < 2e-16 ***
## log(entorno_plazas)         -2.345e-03  2.279e-03  -1.029   0.3038    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4972 on 1127 degrees of freedom
## Multiple R-squared:  0.7628, Adjusted R-squared:  0.7605 
## F-statistic: 329.4 on 11 and 1127 DF,  p-value: < 2.2e-16

Como se explicó, no tiene sentido analizar directamente los parámetros estimados por le modelo lineal sin comprobar previamente la existencia de dependencia espacial. Para ello, es necesario definir una matriz de vecindario que permita incorporar el efecto espacial al análisis. En este caso se define como vecinas a todas aquellas obsevaciones que se encuentran dentro de un radio de 500 metros lineales. Además, se pondera el efecto de cada observación vecina por la distancia, de manera tal que las observaciones más cercanas tengan una influencia mayor que las observaciones más distantes.

library(spdep)

## Cargando paquete requerido: spData

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

library(sp)
cord <- st_coordinates(st_centroid(datos))
d <- dnearneigh(cord, 0, 500)

## Warning in dnearneigh(cord, 0, 500): neighbour object has 298 sub-graphs

dlist <- nbdists(d, coordinates(cord))
idlist <- lapply(dlist, function(x) 1/x)
w <- nb2listw(d, glist=idlist, style="W", zero.policy = TRUE)
w

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 1139 
## Number of nonzero links: 4352 
## Percentage nonzero weights: 0.3354605 
## Average number of links: 3.820896 
## 121 regions with no links:
## 27, 29, 39, 69, 87, 97, 101, 105, 115, 122, 133, 135, 173, 191, 218,
## 219, 224, 226, 228, 229, 232, 249, 260, 263, 279, 292, 298, 302, 306,
## 314, 320, 339, 346, 349, 350, 356, 358, 361, 372, 380, 422, 467, 473,
## 478, 483, 488, 499, 511, 548, 555, 558, 602, 607, 630, 672, 673, 680,
## 682, 683, 701, 705, 713, 730, 731, 749, 751, 778, 779, 780, 784, 792,
## 809, 862, 918, 922, 923, 931, 940, 952, 954, 955, 956, 981, 985, 990,
## 992, 1032, 1045, 1049, 1059, 1068, 1070, 1071, 1072, 1073, 1077, 1080,
## 1081, 1086, 1087, 1089, 1095, 1098, 1103, 1105, 1106, 1108, 1111, 1117,
## 1118, 1120, 1122, 1125, 1126, 1127, 1129, 1131, 1132, 1135, 1136, 1138
## 298 disjoint connected subgraphs
## 
## Weights style: W 
## Weights constants summary:
##      n      nn   S0       S1       S2
## W 1018 1036324 1018 906.7161 4243.495

Una vez construida la matriz \(W\) se procede a realizar las pruebas de hipótesis de multiplicadores de Lagrange robustos para testear si \(\rho=0\) y/o \(\lambda=0\).

lm.RStests(ols, w, test= c("RLMlag","RLMerr"), zero.policy = T)

## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = form, data = datos)
## test weights: w
## 
## adjRSlag = 5.0708, df = 1, p-value = 0.02433
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = form, data = datos)
## test weights: w
## 
## adjRSerr = 629.96, df = 1, p-value < 2.2e-16

Como se puede observar, existe dependencia espacial tanto en la variable dependiente como en el término de error. Por lo tanto, se estima un modelo SAC, que incorpora estos efectos de manera simultánea:

\[Y = \rho W Y + \beta X + \mu\]

\[\mu = \lambda W \mu + \epsilon\]

library(spatialreg)

## Cargando paquete requerido: Matrix

## 
## Adjuntando el paquete: 'spatialreg'

## The following objects are masked from 'package:spdep':
## 
##     get.ClusterOption, get.coresOption, get.mcOption,
##     get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption

sac = sacsarlm(ols, datos, w, zero.policy = TRUE)
summary(sac, Nagelkerke = TRUE)

## 
## Call:sacsarlm(formula = ols, data = datos, listw = w, zero.policy = TRUE)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.1004207 -0.1502080 -0.0086513  0.1322335  2.4949507 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                                Estimate  Std. Error  z value  Pr(>|z|)
## (Intercept)                  8.4157e+00  2.0821e-01  40.4193 < 2.2e-16
## log(dist_bancos)            -2.2793e-01  2.3047e-02  -9.8895 < 2.2e-16
## log(dist_barrio_cerrado)    -2.0360e-02  1.8228e-03 -11.1694 < 2.2e-16
## log(dist_barrios_populares)  2.6146e-04  5.4836e-03   0.0477    0.9620
## log(dist_industrias)        -7.8252e-03  3.8236e-03  -2.0466    0.0407
## log(dist_vias_primarias)    -1.1765e-01  1.2661e-02  -9.2923 < 2.2e-16
## log(dist_vias_secundarias)  -7.8029e-02  1.1442e-02  -6.8196 9.127e-12
## entorno_cuentas              1.3615e-04  1.4246e-05   9.5573 < 2.2e-16
## log(entorno_edificaciones)  -6.4790e-02  1.0547e-02  -6.1429 8.101e-10
## log(entorno_fot)             3.4924e-01  2.4322e-02  14.3588 < 2.2e-16
## entorno_parcelas            -4.1778e-04  5.8737e-05  -7.1128 1.137e-12
## log(entorno_plazas)         -3.8760e-04  1.9862e-03  -0.1951    0.8453
## 
## Rho: 0.017803
## Asymptotic standard error: 0.0097407
##     z-value: 1.8277, p-value: 0.067599
## Lambda: 0.74908
## Asymptotic standard error: 0.01315
##     z-value: 56.966, p-value: < 2.22e-16
## 
## LR test value: 986.4, p-value: < 2.22e-16
## 
## Log likelihood: -321.0514 for sac model
## ML residual variance (sigma squared): 0.07477, (sigma: 0.27344)
## Nagelkerke pseudo-R-squared: 0.90022 
## Number of observations: 1139 
## Number of parameters estimated: 15 
## AIC: 672.1, (AIC for lm: 1654.5)

Como se puede observar, todas las variables resultan estadísticamente significativas, a excepción de la distancia a barrios populares y el porcentaje del entorno cubierto por plazas o parques Sin embargo, la interpretación de los coeficientes no es directa. Esto se debe a que pueden existir efectos cruzados, o “spill-over” entre las diferentes variables independientes y su efecto sobre la variable dependiente. Para poder interpretar correctamente los parámetros estimados es necesario calcular los efectos directos, indirectos y totales, tal cual lo sugiere (Golgher and Voss 2016). Para ello:

impacts(sac, listw = w)

## Impact measures (sac, exact):
##                                    Direct      Indirect         Total
## log(dist_bancos)            -0.2279543528 -4.104282e-03 -0.2320586349
## log(dist_barrio_cerrado)    -0.0203624944 -3.666235e-04 -0.0207291179
## log(dist_barrios_populares)  0.0002614906  4.708098e-06  0.0002661987
## log(dist_industrias)        -0.0078261456 -1.409085e-04 -0.0079670541
## log(dist_vias_primarias)    -0.1176635385 -2.118513e-03 -0.1197820520
## log(dist_vias_secundarias)  -0.0780377901 -1.405058e-03 -0.0794428482
## entorno_cuentas              0.0001361674  2.451673e-06  0.0001386191
## log(entorno_edificaciones)  -0.0647977166 -1.166673e-03 -0.0659643893
## log(entorno_fot)             0.3492810038  6.288749e-03  0.3555697532
## entorno_parcelas            -0.0004178321 -7.523001e-06 -0.0004253551
## log(entorno_plazas)         -0.0003876468 -6.979519e-06 -0.0003946263

Por lo tanto, por ejemplo, se observa que un aumento de un 1% del factor de ocupación de la tierra (entorno_fot) se traduce en un aumento del 0.34% en el precio de la tierra urbana. De igual manera, un aumento de un 1% en la distancia a una via primaria (dist_vias_primarias) implica una reducción del precio de la tierra del 0.12%. También se puede observar que el precio de la tierra es mayor mientras la menor sea la distancia a barrios cerrados (dist_barrio_cerrado). El resto de los coeficientes se puede interpretar de la misma manera. Estos impactos son de utilidad para la valuación de proyectos de infraestructura urbana, en donde las externalidades generadas suelen ser capturadas sin contrubución por parte de algunos agentes económicos.

Introducción a la aplicación de técnicas de aprendizaje computacional a fenómenos territoriales

En esta sección se analizará brevemente la librería de R “valuate” (Carranza 2025). Esta librería provee una serie de funciones orientadas a facilitar el proceso de valuación de inmuebles, siguiendo el desarrollo de (carranza_eguino?) (en prensa). Sin embargo, puede ser también de utilidad para la predicción de la distribución espacial de cualquier tipo de fenómeno geográfico.

El flujo de trabajo propuesto apunta a la generación de herramientas para facilitar las siguientes acciones:

• Creación de variables territoriales que informan sobre la distancia a diferentes hitos o las características de determinados entornos.
• Homogeneización de valores para descontar los efectos de características intrínsecas de las observaciones. muestrales. Por ejemplo, en el caso de la valuación de la tierra, se descuentan los efectos sobre el precio por metro cuadrado de la forma de las parcelas, su ubicación en la cuadra o el tamaño del lote.
• Diferentes alternativas para la identificación de observaciones muestrales espacialmente atípicas.
• Entrenamiento de diferentes modelos de aprendizaje computacional y estimación del nivel de error esperado fuera de la muestra.
• Simulación de escenarios que permiten conocer la magnitud y extensión geográfica de diferentes externalidades generadas por modificaciones en las variables independientes sobre el precio de los inmuebles.

Instalación

Para instalar la versión de desarrollo de esta librería en la consola de R se debe ejecutar el siguiente código:

remotes::install_github("carranzajuanp/valuate")

Contenidos

La función calcular_dist() permite generar un ráster que informa la distancia de cada píxel hacia el hito más cercano.

Por ejemplo, si se desea conocer la distancia a la vía primaria más cercana en la ciudad de La Plata (Argentina):

# Cargamos la librería al entorno de trabajo
library(valuate)

# Definimos el área de estudio
area_de_estudio <- nominatimlite::geo_lite_sf(address = "La Plata, Argentina", points_only = F)
bbox = sf::st_transform(area_de_estudio, 4326) |> sf::st_bbox(bbox)

# Descargamos las vias principales
vias_prim <- osmdata::opq(bbox) |>
  osmdata::add_osm_feature(key = "highway", value = "secondary") |>
  osmdata::osmdata_sf()
vias_prim <- vias_prim$osm_lines

# Calculamos un ráster con las distancias a la vía principal más cercana
calcular_dist(area = area_de_estudio,
              objeto = vias_prim,
              dim = 50,
              nombre = "vias")

## El proceso de calculo ha finalizado, y el raster resultante fue guardado en el directorio de trabajo con el nombre vias.tif. En el environment se ha creado el raster 'vias'.

La función calcular_entorno() resume las características de un entorno definido por el usuario, a partir de capas geográficas vectoriales de polígonos o de puntos. Cuando se le provee un objeto espacial “sf” con geometría “POLYGON” o “MULTIPOLYGON” la función permite calcular variables independientes que resumen las características del vecindario, indicando la proporción del espacio que se encuentra cubierto por los polígonos de interés. Por ejemplo, permite calcular el porcentaje del vecindario que se encuentra cubierto por plazas o parques. Cuando se le provee un objeto espacial “sf” con geometría “POINT” la función permite calcular variables independientes que resumen las características del vecindario, indicando la cantidad de puntos de interés que se encuentran dentro del entorno definido por el usuario. Por ejemplo, permite calcular la cantidad de negocios u oficinas que se encuentran dentro de un radio de “x” metros de cada píxel.

Por ejemplo, si se desea conocer información sobre los espacios verdes en el vecindario en la ciudad de La Plata (Argentina):

# Definimos el área de estudio
area_de_estudio <- nominatimlite::geo_lite_sf(address = "La Plata, Argentina", points_only = F)
bbox = sf::st_transform(area_de_estudio, 4326) |> sf::st_bbox(bbox)

# Descargamos los parques que se encuentran dentro del área de estudio
plazas <- osmdata::opq(bbox) |>
  osmdata::add_osm_feature(key = "leisure", value = "park") |>
  osmdata::osmdata_sf()
plazas <- plazas$osm_polygons

# Calculamos un ráster que resume, para cada píxel las características del entorno
calcular_entorno(area = area_de_estudio,
                 ext = 500,
                 objeto = plazas,
                 dim = 50,
                 nombre = "plazas_entorno")

## Warning in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy, :
## GDAL Error 1: PROJ: proj_as_wkt: DatumEnsemble can only be exported to
## WKT2:2019

## [1] "El proceso de calculo ha finalizado, y el raster resultante fue guardado en el directorio de trabajo con el nombre plazas_entorno.tif. En el environment se ha creado el raster 'plazas_entorno'."

La función calcular_raster() devuelve un ráster que resume el cálculo realizado a partir de otro ráster, incluyendo fuentes de datos con diferente proyección y resolución espacial. La función permite readecuar las resoluciones y proyecciones de diferentes rásters a las utilizadas en el proyecto. Además, el paŕametro “entorno” permite realizar un cálculo del promedio de la variable bajo análisis en el vecindario definido por el usuario. Por ejemplo, si se desea readecuar a los criterios de un proyecto un ráster que informa sobre el porcentaje de edificaciones detectadas en un píxel a partir de una imagen satelital en la ciudad de Medellín (Colombia):

original = terra::rast("edificaciones.tif")
area_de_estudio = sf::st_read("area.gpkg")

## Reading layer `area' from data source 
##   `/home/juan/Documentos/GitHub/spatial_data_science/area.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -75.66658 ymin: 6.066645 xmax: -75.49011 ymax: 6.367076
## Geodetic CRS:  WGS 84

calcular_raster(original, area_de_estudio, dim = 50, entorno = 100, "prueba_raster")

## |---------|---------|---------|---------|=========================================                                          

## El proceso de calculo ha finalizado, y el raster resultante fue guardado en el directorio de trabajo con el nombre prueba_raster.tif. En el environment se ha creado el raster 'prueba_raster'.

La función entrenar_modelo() permite entrenar modelos y cuantificar su error en la predicción del valor de los inmuebles (o cualquier otra variable territorial que se esté analizando). Esta función implementa un flujo completo de entrenamiento, evaluación y ajuste de un modelo supervisado para la estimación de una variable dependiente a partir de un conjunto de variables independientes en formato ráster. Específicamente, está diseñada para trabajar con datos espaciales almacenados en un dataframe con geometría de tipo ‘POINT’, usando herramientas de las librerías ‘sf’ (Pebesma et al. 2023) y ‘terra’ (Hijmans 2024). El proceso se desarrolla en los siguientes pasos:

• Preprocesamiento de datos: Los rásters de las variables independientes se combinan con el dataframe espacial mediante la extracción de valores correspondientes a los puntos de la geometría. Los datos se depuran eliminando filas con valores faltantes.

• Configuración del modelo: Se define un modelo a través de la librería ‘caret’ (Kuhn et al. 2023) , con soporte para una amplia variedad de algoritmos. El modelo predeterminado es Quantile Regression Forest (QRF), aunque se pueden especificar otros métodos soportados por ‘caret’.

• Entrenamiento iterativo: La función utiliza validación cruzada de 10 particiones para evaluar el error medio absoluto porcentual (MAPE) fuera de la muestra. Si el MAPE promedio supera el umbral especificado, se eliminan observaciones cuyo MAPE individual sea igual o superior al 40%, y el modelo se entrena nuevamente con los datos restantes. Este proceso se repite hasta que el MAPE promedio cumpla con el umbral definido.

• Generación de predicción espacial: Una vez finalizado el entrenamiento, el modelo entrenado se utiliza para interpolar las predicciones sobre los rásters de entrada, generando un nuevo ráster de salida que se guarda en el directorio de trabajo con el nombre ‘vut.tif’.

• Salida de datos: La función genera dos conjuntos de datos en el entorno global: datos_utilizados (datos que cumplieron con el umbral de error) y datos_eliminados (datos descartados por errores altos).

La función está optimizada para aplicaciones que requieren análisis espaciales precisos y está diseñada para integrarse en flujos de trabajo geoespaciales complejos. Además, emplea paralelización para acelerar el proceso de entrenamiento en máquinas con múltiples núcleos. El diseño de esta función está orientado a contextos donde es crucial garantizar una alta precisión en la estimación, especialmente en aplicaciones que involucran modelización espacial o análisis urbano. Además, facilita la integración con flujos de trabajo basados en datos geoespaciales mediante las librerías ‘sf’ y ‘raster’.

# A continuación se presenta un ejemplo para la Ciudad de Medellín (Colombia), con sólo 500 datos.
load("datos_ml.rda")

entrenar_modelo(df = dat,
                dependiente = "vut",
                independientes = c("dist_basura.tif",
                                   "dist_industria.tif",
                                   "dist_plazas_parques.tif",
                                   "dist_tren.tif",
                                   "dist_vias_autopista.tif",
                                   "dist_vias_prim.tif",
                                   "dist_vias_sec.tif",
                                   "entorno_altura.tif",
                                   "entorno_plazas_parques.tif",
                                   "entorno_comercios.tif",
                                   "entorno_edificaciones.tif",
                                   "entorno_hoteles.tif"),
                modelo = "qrf",
                umbral = 0.3)

## Reading layer `area' from data source 
##   `/home/juan/Documentos/GitHub/spatial_data_science/area.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -75.66658 ymin: 6.066645 xmax: -75.49011 ymax: 6.367076
## Geodetic CRS:  WGS 84

## Number of pixels is above 5e+05.Only about 5e+05 pixels will be shown.
## You can increase the value of `maxpixels` to 1651724 to avoid this.

Finalmente, la función simular_escenario() permite, a partir de la predicción de la distribución espacial de una variable, estimar escenarios a partir de la simulación de alteraciones en alguna (o algunas) de las variables independientes. La función devuelve un ráster que informa sobre el impacto en la variable dependiente de las simulaciones realizadas en las variables independientes. Permite estimar la magnitud y el alcance territorial de las modificaciones simuladas sobre la variable de estudio.

Referencias

Anselin, Luc. 1988. Spatial Econometrics: Methods and Models. Vol. 4. Studies in Operational Regional Science. Dordrecht: Springer Netherlands. https://doi.org/10.1007/978-94-015-7799-1.

———. 1995. “Local Indicators of Spatial Association LISA.” Geographical Analysis 27 (2): 93–115. https://doi.org/10.1111/j.1538-4632.1995.tb00338.x.

Appelhans, Tim, Florian Detsch, Christoph Reudenbach, and Stefan Woellauer. 2023. “Mapview: Interactive Viewing of Spatial Data in r.” https://CRAN.R-project.org/package=mapview.

Carranza, Juan Pablo. 2025. Valuate: Caja de Herramientas Para La Realización de Estudios Territoriales. https://github.com/carranzajuanp/valuate.

Golgher, André Braz, and Paul R. Voss. 2016. “How to Interpret the Coefficients of Spatial Models: Spillovers, Direct and Indirect Effects.” Spatial Demography 4 (3): 175–205. https://doi.org/10.1007/s40980-015-0016-y.

Hernangómez, Diego. 2024. “Nominatimlite: Interface with Nominatim API Service.” https://doi.org/10.32614/CRAN.package.nominatimlite.

Hijmans, Robert J. 2023. “Raster: Geographic Data Analysis and Modeling.” https://CRAN.R-project.org/package=raster.

———. 2024. “Terra: Spatial Data Analysis.” https://CRAN.R-project.org/package=terra.

Hijmans, Robert J., Márcia Barbosa, Aniruddha Ghosh, and Alex Mandel. 2024. “Geodata: Download Geographic Data.” https://CRAN.R-project.org/package=geodata.

Kuhn, Max, Jed Wing, Steve Weston, Andre Williams, Chis Keefer, Allan Engelhardt, and Tony Cooper. 2023. Caret: Classification and Regression Training. https://cran.r-project.org/web/packages/caret/index.html.

Manski, Charles. 1993. “Identification of Endogenous Social Effects: The Reflection Problem.” The Review of Economic Studies 60 (3): 531–42. https://EconPapers.repec.org/RePEc:oup:restud:v:60:y:1993:i:3:p:531-542.

OpenStreetMap, contributors. 2023. “Planet Dump Retrieved from Https://Planet.osm.org.” https://www.openstreetmap.org.

Padgham, Mark, Bob Rudis, Robin Lovelace, and Maëlle Salmon. 2017. “Osmdata.” https://doi.org/10.21105/joss.00305.

Pebesma, Edzer, Roger Bivand, Etienne Racine, Michael Sumner, Ian Cook, and Tim Keitt. 2023. Sf: Simple Features for r. https://cran.r-project.org/web/packages/sf/index.html.

Wickham, Hadley, Romain François, Lionel Henry, Kirill Müller, and Davis Vaughan. 2023. “Dplyr: A Grammar of Data Manipulation.” https://CRAN.R-project.org/package=dplyr.